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The stability of the ideal magnetohydrodynamic (MHD) interchange mode at marginal conditions is studied. 
A sufficiently strong constant magnetic field component transverse to the direction of mode symmetry provides 
the marginality conditions. A systematic perturbation analysis in the smallness parameter, \bij B^l 2 , is 
carried out, where B c is the critical transverse magnetic field for the zero-frequency ideal mode, and 62 is the 
deviation from B c . The calculation is carried out to third order including nonlinear terms. It is shown that 
the system is nonlinearly unstable in the short wavelength limit, i.e., a large enough perturbation results in 
instability even if b 2 /B c > (linearly stable). The normalized amplitude for instability is shown to scale as 
\i>2/ B c \ x / 2 . A nonlinear, compressible, MHD simulation is done to check the analytic result. Good agreement 
is found, including the critical amplitude scaling. 



I. INTRODUCTION 

It is well known that the magnetohydrodynamic 
(MHD) magnetized plasma interchange instability can 
be stabilized by a transverse magnetic field. For a given 
wavenumber, allowing a magnetic field component in the 
direction of the wavenumber introduces Alfvenic stabi- 
lizing tension such that beyond a critical transverse field 
(transverse to the direction of mode symmetry) that 
wavenumber is linearly stableP The nonlinear evolution 
of the magnetized plasma interchange instability is less 
well understood. In particular, the state of the system for 
when the transverse B-field is marginally subcritical (or, 
equivalently, the plasma beta is slightly above critical) is 
an important question for magnetized fusion energy ap- 
plications: does the mode saturate at low amplitude and 
how does the marginal convection and resulting transport 
scale with deviation from marginality? The question is 
an important consideration for stellarators, for example, 
since these fusion devices are engineered for very high 
precision magnetic fields and one of the precision con- 
straints arises from ideal MHD linear stability resultsP If 
the nonlinear consequence of a slightly subcritical B-field 
were better understood, it may be possible to optimize 
over the MHD design constraints. It was also recently 
shown that the linear growth rate of ideal interchanges 
in a reversed-field pinch for a slightly subcritical B-ficld is 
weaker than expected and may be overcome by nonlinear 
effect sP 

The interchange mode is a pressure-driven mode that 
is characterized by the interchange of magnetic flux tubes 
so that the overall free energy of the system is loweredP 
The instability occurs when the equilibrium has a density 
gradient unfavorable to the direction of a "gravitational" 
force. In systems with curvature, this force comes from 
a centrifugal force generated by thermal motion in field 
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curvature. The mode can be stabilized by introducing a 
strong enough field, transverse to the flutes, that prevents 
the flux tubes from being able to freely interchange. The 
strength of the stabilizing field can be determined using 
linear theory and will depend on the steepness of the 
density gradient and the magnitude of the gravitational 
force. 

There have been a few studies done on nonlinear 
growth of interchange instabilities at marginal stability 
in tokamaksPESl Although a Lagrangian approach has 
been attemptedp-^the general approach is to expand the 
equation of motions of the unstable mode about marginal 
stability and thus the nonlinear terms in the system can 
be evaluated!^ We can determine the overall stability 
of the system by comparing the behaviour of the non- 
linear effects to the linear driving term. In Ref. [5] the 
author found that, for the profiles investigated, the non- 
linear effects were stabilizing. Similarly, in Ref. ;7: the au- 
thor showed nonlinear saturation at marginal stability. 
Both authors considered a system with a sheared mag- 
netic field. In studying the line-tied g mode, the authors 
in Ref. 8s showed that near the marginally stable point the 
system was nonlinearly unstable. However, Refs. [9] and 
1101 showed that the nonlinear growth transitions through 
an initial regime where the nonlinear growth dominates 
the linear response, as shown in Ref. but a secondary 
regime takes over when the amplitude is sufficiently large 
and so the mode amplitude remains bounded. 

We simplify our system to a slab geometry where we 
use an effective gravitational field, g, to model centrifu- 
gal force due to field line curvature^ and we assume a 
constant transverse field. This reduces the complexity of 
the system so that the focus of the analysis can be on 
how nonlinear terms get introduced into the equations of 
motion. The idealized system is described in Sec. |ll| along 
with the derivation of nonlinear time evolution equation. 
The goal is to have a simpler methodology in a simple 
system that can be generalized into more complicated 
systems, e.g. sheared fielcP^I, ballooning * 13 * 14 !, etc. In 
Sec. |III| we verify our result using a dissipative numer- 
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ical simulation. The results are summarized in Sec. IIV 



II. THEORY 

Consider a slab system with constant gravity g = — gx 
and very strong magnetic field in the z-direction such 
that B± <C B z and Vaz is much larger than u, the typical 
flow in the system. This system is incompressible and 
can be described by the two-dimensional MHD reduced 
equations^ given by, 



d t p + u • Vj_p = 0, 



z • V ± x (p-u) 



B • VxViv + gdyP, 



dtip — B ■ Vx< 



0. 



(1) 



(2) 



(3) 



where u = i x V±<f> and Bj^ = z x and we have 

defined, in the usual way, 



d 
dt 



d_ 
Ft 



u- Vx- 



Variations in z are suppressed since the fastest inter- 
change has d/dz = 0. The nonlinear system of equations 
Q-Q can be solved for the variables p, the density, and 
<j> and ip, the flow and magnetic streamfunctions, respec- 
tively. 

We consider a static equilibrium with a density gradi- 
ent that's unstable to interchange and a constant, trans- 
verse magnetic field. More explicitly, we have 



p' (x) > 0, B = B y, 



0. 



(4) 



where henceforth the primes denote differentiation with 
respect to x. We also add the assumption that p' Q —> 
at the boundaries and has even parity. 

Small perturbations about this equilibrium yield the 
WKB dispersion relation 



k 2 Vl y 



7 S 2 



(5) 



where 7 g = Igp'/p] 1 ^ 2 is the Rayleigh- Taylor growth 
rate and k is the wavenumber in the y direction. 
The Rayleigh- Taylor growth rate represents the effective 
"gravitational" acceleration and is the driving force in 
an interchange instability. With a transverse magnetic 
field, field line bending results in Alfvenic restoring forces 
with frequency kVAy In this paper, we consider the dy- 
namics of the magnetized interchange mode when the 
magnetic field strength is strong enough to just stabilize 
interchanges, i.e., the system is near marginal stability. 
In particular, for a given k, suppose w 2 > everywhere 
in x except for a single small region where it is very close 
to zero, positive or negative. In that case, weakly grow- 
ing perturbations are possible in the vicinity of where 
k V? — 7g is close to zero. The time rate of change 



of the perturbations will be very small compared to the 
local 7 9 . Thus, we order 



e< 1. 



(G) 



This implies that any deviations in Bq away from criti- 
cality must be small. In particular, if 



B = B C + b 2 



(7) 



then, according to 62/^ must be of 0(e 2 ). 

We allow small perturbations about this marginal 
point such that the amplitude of the magnetic perturba- 
tion, A, while small, is large enough that the nonlinear 
magnetic tension forces can influence the growth time. 
This results in the optimal ordering 



(8) 



We represent the perturbation by expanding ip and (f> in 
a series 



(9) 



(10) 



where the order in e is denoted by the subscript. The 
continuity equation ([!]) can be satisfied by letting p — 
p(ip) and using With this change of variable we can 
expand p in terms of Sip — 4> — to get 



pW = p° + jt-W 

Bq 



Substituting the expansions ©-gl) into Q and g we 
can solve for the nonlinear evolution of the perturbations 
order by order. 



A. First order equations 

Matching terms to lowest, non- vanishing order, we ob- 
tain from ([2]) and ^ the equations 



= B c d v V±tpi + gdypx, 



where (11) gives 



B c dy<f>! = 0, 



Pi = Po^i/Bc. 



Substituting p\ into (12) the equation becomes 



where we have defined the operator 

£(/) = (Vi + ^ P ' )dyf. 



(12) 
(13) 

(14) 

(15) 
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Writing tpx as 

1> 1 (x,y,t) = A(t)( 1 (x)caB(ky), (16) 
we obtain the eigenvalue equation 

C('(x)-fc 2 Ci(x) + ±fai(x) = (17) 

that can be solved to get an eigenvalue for B c . The 
boundary condition for po implies that £i decays expo- 
nentially close to the boundary. 

In writing ( |16[ ) we assumed a cosine perturbation in 
the density which implies ipi ~ cos(fcy) from (111. If 



we also assume that this initial perturbation results in a 
pure mode for the lowest order flow then 



<j>i(x,y,t) = 



(18) 



is the solution to (131 



To the lowest order we have found that given the mode 
of the density perturbation, fc, and the equilibrium den- 
sity gradient profile, p' {x), then the marginally stable 
field strength B c can be solved for using (17 1. This re- 



sult is consistent with the prediction from linear theory 
for the existence of the marginally stable value. 



B. Second order equations 

In order to solve for the time evolution of ipi, it is 
necessary to proceed to higher order in the expansion. 
We now match 0(e 2 ) terms in equations ^ and ^ to 
get 

B c d y Vi^ 2 + gd yP2 + Bi ■ V ± Vi^i = 0, (19) 



d t ipi = B c d y (j)2, 



where (111 to the same order gives 



(20) 



(21) 



Using ipi from ( 16 1, 4> 2 can be solved for in ( 20 1 to obtain 
<h{x,y,t) = 1 dA jt) ft (g) sin(fc2/) + fc( Xj t)i (22) 



kB c dt 

where 4> 2 (x,t) is a constant of integration. 

Substituting (21) into (19) results in an equation for 

■02, 



£(02) 



m 



Po d v($i)> 



(23) 



where we have used (151 to simplify the Laplacian. This 



has a solution of the form 

ih{x, y, t) = A(t) 2 ( 2 (x) cos(2fcy) + $ 2 (x, t), (24) 



where ip2 (%, t) is the homogeneous solution to (23 1. Sub- 
stituting p4| into ( 23 ) , we find £2 (%) satisfies 



(^ x )-4k\ 2 (x) + ^p' ( 2 (x) 



2B? 



PoCi(^) 2 - (25) 



To fully analyze the stability of our system we still have 
to resolve the time-evolution of ipi . It is also important 
to solve for ip 2 and 4>2 to make sure that those terms are 
well-behaved. 



C. Third order equations 

As was done previously in lower orders, we match 
terms of 0(e 3 ) in ^ and pfc. The resulting higher order 
equations are 

<9t(poVi(fc + p' cj>' 2 ) = BcfyViVs + b 2 d y Vl^! 
+ gd yP3 +B 1 ■ V x ViV2+B 2 - VxViVi (26) 



d t <if) 2 = B c d. 



■ Bi • Vj_< 



along with 



P3 



Po , P0 b 2 , 

B~ r ^ - ~bJ^ 



PL 

Bl 



1 n'" 

6 b^ 1 



(27) 



(28) 



from (111. 



Integrating ( 26 ) over one period in y we find that <p 2 



is not driven by tpx so we can set 

fa(x,t) =0 



(29) 



without loss of generality. No zonal flows are generated 
in the system when creating a periodic perturbation in 
the density. However, averaging (27 1 over y we find that 



zonal fields are generated according to 

U^t) = ~A{tfQ l {x)C[{x). 



(30) 



For a given k and po the system is now solved up to 
second order with the exception of the time-evolution 
A(t). The variables tpx, (f>i, ip 2 , and (f> 2 are defined by 



(161, (18), (24), and (22), respectively. We can solve 
for £1 using (17) and then for ( 2 using (25|. The y- 



independent terms 4> 2 and ip 2 are given by ( 29| and ( 30 1 . 

To solve for A(t) we need to simplify ( |26[ ) by making 
use of ( 15 1 , ( 20 ) , and ( 28 ) . After some algebra ( 26 ) takes 



the form 
1 



k 2 B d ?(# 2 -PoPo d y'4>i ~ Po d M 



2-^6 2 Po^i+^[0i>2] 



(31) 



where exact details of the functional T is suppressed here 
for clarity but is shown in Appendix A. We can extract 



a time-evolution equation by substituting (16), (24), and 
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( 30 1 into the above equation and applying the operator 
Jdx (i(x)fd(cos(ky)) evaluated over all space. This op- 
eration will annihilate the ip3 term and any higher order 
harmonics. 

After simplification (see Appendix A), we arrive at the 
equation for A(t) 

^4"(Pop[,Ci 2 >^^W = -2b 2 (p' ( 2 )A(t) 



+ ( <Pod 2 C 2 ) - ^U,C?(C?)"> 



IB, 
"4 9 



(Ci 2 (Ci 2 )"") A{t)\ 



(32) 



The time-evolution equation ( 33 ) closes the system and 



we can fully determine the first and second order pertur- 
bations, ipi and ip2 defined by (16) and (24), for a given 
k, po, and b 2 . This is achieved by first solving the eigen- 



value problem ( 36 ) and using the solution for Z\ and B c 
to solve for Z 2 using (37), and finally determining the 



coefficients (34), (35) and solving for A(t) in (33). 

The coefficient c\ is a positive number for p' > 0, and 
so the linear stability of the system is determined by the 
sign of b 2 • This result agrees with the linear theory. How- 
ever, the overall nonlinear stability of the system is going 
to be determined largely from the sign of C3 compared to 
the sign of b 2 . 



where the angled brackets are defined as 

(J) = ^Jdxf(x) 

with L~ l = p'q/po evaluated at x = 0. We can simplify 
this further by letting 

x -> xL p , 
p (x) -> po(0)p(x), 

-> Zi(x), 
( 2 (x) -> Z 2 ( X )/(B c L p ), 

in order to introduce dimensionless versions of the vari- 
ables x and po, and have A with dimensions of ip. Ap- 
plying this normalization to ( 32 ) we get 



1 



k 2 V 2 Ac dt 2 



A(t) = 



-2^-c x A{t) 



C3 



B 2 C L 2 



where V\ c = 



B 2 c /p {0) and 



c\ = 



(p'zr, 



A{t) 3 (33) 



(34) 



C3 



(p"ZlZ 2 ) 

WzJ) 



1 {p'ztjzt)") 

4 {pp'Z 2 ) 

1 V 2 C (Zf(Zf)"") 
4 gL p {pp'Z 2 ) 



(35) 



where the primes and brackets now denote derivatives 



and integrals in \- Using the same normalization on ( 17) 
and (25) we get the following equations, 



z';-k 2 L 2 p z 1 + 9 ^fp' Zl = ^ 



Z'' - 4k 2 L 2 Z 2 



gL P , „ 1 gLp „ 2 



n" 7 l 



-p Z 2 = -- 



(36) 



(37) 



for the dimension- free Z\ and Z 2 . 



D. Short wavelength limit 

We can analytically solve ( [36] ) for the case kL p 3> 1 in 
which regime the cells are elongated in ^-direction but 
still shorter than the scale of the gradient, i.e., 

kL p > > 1. 

With this scaling we can approximate p'(x) to be 



Ax) 



2 



(38) 



Assuming that gL p /Vj : 



k L p , then from scaling ar- 



guments we find that ( 36 1 has the familiar form of a 



quantum harmonic oscillator. This has the well-known 
solution 



z i (x) = z i exp 



kL 



p 2 
X 



k 2 L 2 = 



gLp 



1 - 



2V2 



1 1 

V^kLp 



(39) 



(40) 



for the ground state. 
kLp ^> 1 and the solution for the 



This solution is correct only for 
energy" adds a small 
correction to the initial assumption. Using the same scal- 
ing, to lowest order, (37) has the solution 



z 2 (x) = -\xZi{xY l - 



(41) 



The time-evolution equation ( 33 ) can be simplified in 
the kLp 1 limit by substituting the solutions (39 (-(41 ) 
in the coefficients ( 34 ) and ( 35 1 . After simplification wc 



arrive at the following values for the coefficients 



ci 



1. 



1 



C3 



-kL n 



(42) 



where we only kept the largest terms and have assumed 
that Z x = 1. 

The above result implies that even if b 2 > 0, if the 
initial amplitude A = A(0) is such that 



An 



BrLr 



> 4, 



1 



Br kL„ 



(43) 
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then the system will be nonlinearly unstable and the am- 
plitude will increase without bound. Furthermore, for 
62 < the instability grows faster than predicted from 
linear theory and any small perturbation will continue to 
grow larger without saturation. 

With the solution for the eigenmode we can check the 
ratio between the spatial scale of the perturbation, char- 
acterized by the displacement in the x-direction £ x , and 
the width of the eigenmode 



A 



(44) 



given by ( 39 ) . The displacement is related to the velocity 
such that dt£,x ~ u x , and from (20) we get that dtA ~ 
B c u X 2, which implies that A ~ B c ^ x . Substituting for A 



using (43 1 gives us a scale for the displacement, 



B r k 



which yields 



(45) 



(46) 
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for the ratio of the two scale lengths. As should be ex- 
pected, the spatial size of the amplitude required to be 
nonlinearly unstable is much smaller than the width of 
the eigenmode. 



FIG. 1. The equilibrium profiles for the background field 
B z , the density p, and the magnetic streamfunction ip along 
with the difference from constant field. 



III. NUMERICAL SIMULATION 

To confirm this result, we used a two-dimensional code 
that solves the fully compressional equations (see Ap- 
pendix B). The variables p, pu, ip, and B z are solved 
numerically and stepped in time. We set B z ^> \B±] so 
the equations are effectively reduced. The code is dissi- 
pative so we introduced source terms in the density in 
order to maintain a steady state profile suitable for our 
model. The sourcing, although weak, results in a profile 
for By(x). To compare with analytic theory, we wish to 
keep By approximately constant. Thus, we allowed B z 
to resistively relax at a somewhat slower rate than B y in 
the equilibrium. 

The system is normalized so that initially Vaz — 1 and 
L x = 1, where L x is the height of the box. We used 
hardwall, free-slip boundary conditions for the top and 
bottom walls and periodic boundary conditions for the 
sides. The periodic boundary conditions discretize the 
system so that the only wavenumbers allowed are integer 
multiples of 2ir/L y , where L y is the width of the box. 
From ^ we know that the lower modes are the most 
unstable, so to study the case with kL p ^> 1, i.e. short 
wavelength, we selected L y such that the minimum value 
for kL p satisfies this condition. By choosing k = 2n/L y 
we can satisfy the marginality condition by adjusting B 
and/or g such that kVA y ~ 7 9 for the minimum mode. 



We set L y — 0.5, and from the density profile we have 
L p = Pq/p'q ~ 0.4 and so we satisfy the condition 

kLp 5.03 > 1 

which is necessary to compare with the analytical re- 
sult from Section |II D| We attempted to run tests with 
a larger value of k by decreasing L y , but the code was 
numerically unstable for smaller box widths. 

To generate the equilibrium we initialize i/i to Bqx and 
let the system reach an equilibrium which is steady state. 
The density source term results in a weak flow in the in- 
direction. This flow scales with the diffusion, so a mini- 
mal, numerically-stable value for the diffusion is chosen 
to minimize its effect. The equilibrium profiles for the 
density and the background field generated are shown in 
Fig. [I] It is important to note that the equilibrium profile 
for the density does not have p' — > at the boundaries. 
The boundary conditions imply that p' — > —gpo at the 
wall. 

After the equilibrium is made, a density perturbation 
is introduced with p(x,y) — ciocos(ky). From ( |14[ ) we 
can relate the density perturbation amplitude, a(t), to 
the perturbation amplitude of tp, i.e. a = PqA/B c . In 
Fig. [2] we show the resulting unstable eigenmode devel- 
oping for the density. For tests done with Bq far away 
from marginality, i.e. I&2/-EM ~ 50%, there was excellent 
agreement for the growth rate/frequency in the simula- 
tion with ([5]). The theory predicts that there will be 
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FIG. 2. The linear growth of an unstable localized mode 
cut at y ~ 0.26 and for t < 60ta- Time traces separated by 
t « 6ta are shown. 



nonlinear coupling to the mode with wavenumber 2k, so 
it is important that this mode and higher modes are al- 
lowed. Since the diffusivity is weak, it is ensured that 
this is the case. 

Since we can adjust both Bq and g to achieve marginal 
stability, we decided to fix the value of g at 0.15, and ad- 
just B . With this value of g we expect that B c « 0.05 
based on ( 40 1 . However, we found that an equilibrium 
with Bq — 0.05 is stable to perturbations as large as 
10 _1 in the simulation. We decreased the strength 



a 

of the transverse field until it became unstable to pertur- 
bations with ao = 10 -4 . This value was at Bo « 0.0438 
and we took this to be the critical value of the trans- 
verse field for the numerical simulation. Since the criti- 
cal amplitude scales like the square root of the deviation 
from marginality, we are limited to perturbations only 
as small as 10~ 4 otherwise smaller perturbations would 
have meant having deviations that are close to the limits 
of our computational power. 

We created multiple equilibria with different transverse 
field strength within 10% of the numerical critical field 
strength. These equilibria were then perturbed with ao 
of different orders of magnitude. The result of the test 
is shown in Fig. [3] where circles and crosses mark stable 
and unstable points, respectively, and the solid line is 
for a = 4p' Q (b 2 /B c )(L p /k), from our theory, using the 
parameters from the numerical simulation. The slope of 
the theory line seems consistent with the numerical data, 
however, the theory requires larger ao for nonlinear in- 
stability. This inconsistency could be due to the diffusion 
in the code and, in particular, the resistivity may allow 
for slippage in the magnetic field lines which can shift the 
stability boundary at marginal stability. We can calcu- 
late the scale size of this shift based on the values used 
in the simulation (see Appendix B), 



FIG. 3. Result of stability test for a range of deviations from 
B c and magnitude of perturbation, ao. Stable and unstable 
results are denoted by a circle or a cross, respectively. The 
solid line is the theoretical boundary. 



This implies that there could be a shift in B c of or- 
der \Jb2/B c . At marginal stability, even small diffusion 
can cause significant shifts in stable-unstable boundaries. 
However, this implies a shift in B c ; it is harder to explain 
why resistivity results in a nonlinear instability at large 
amplitude of perturbation. It is possible that diffusive 
effects may affect the critical amplitude for nonlinear in- 
stability, but the existence of a nonlinear instability phe- 
nomenon is harder to explain as a diffusive effect. 

In addition to checking the perturbations for a growing 
linear mode, we also check the time trace of the ampli- 
tude for nonlinear effects. In Fig. [4] we show a time trace 
of the amplitude of p, a(<), for the same -Bo but differ- 
ent ao- We can see that the behaviours are different for 
the two cases. In the unstable case, Fig. [4^,, the density 
perturbations become very large quickly and eventually 
dissipate after it hits the boundaries (t < lOOr^). The 
time trace of p' shows that the density profile flattens 
out (// — > 0) after reaching a peak. So, even though our 
analysis in Sec. [Tl] is only valid as long as A < e we can see 
from the trace that it continues beyond this limit until 
the profile collapses. The stable case, Fig. [4}}, has an ini- 
tial growth eventually hitting a peak and then has stable 
oscillations. Even though the amplitude increases some, 
it is still small and the density profile holds. This can be 
seen from the fact that p' is staying constant the entire 
time. We can see in Fig. [5] that as we increase b2/B c 
further from marginality, this initial growth decreases in 
magnitude. It also develops faster and has more noise 
that is indicative of a transient oscillatory mode. 



IV. SUMMARY AND CONCLUSIONS 



(47) 



In this paper, we studied the nonlinear behaviour 
of a marginally stable interchange system. We used 
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FIG. 4. Time trace of the amplitude of density perturbations 
p (solid line) and x derivative of the density p' (dashed line) 
for b 2 /B c w 0.04% with (a) a = 10" 2 and (b) a = 10" 4 . 



the reduced equations to find an analytic solution near 
marginality given a density profile, po{x), deviation from 
marginality, 62, and wavenumber of perturbation, fc, of 
the B-field. The result is a nonlinear differential equa- 
tion for the amplitude of the density perturbations as a 
function of time. The threshold for nonlinear instability 
is dependent on the above quantities, along with g. The 
principal finding of this paper is that marginally stable 
interchange modes in a magnetized plasma can be non- 
linearly unstable for large enough initial perturbations. 
We arrived at this result from a systematic asymptotic 
expansion about marginality in the smallness parameter, 
fa/Be] 1 / 2 , carried out to third order. The first order so- 
lution can be found using the linear eigenvalue problem. 
This solution is then used as a source for the second or- 
der problem. The third order analysis yields the equation 
for the time-dependence of the perturbation. We found 
that the stability of the solution can be determined by 
calculating the coefficient of the nonlinear term in the 
differential equation. This is a nontrivial task for a gen- 
eral perturbation, but we could analytically solve this in 
the short wavelength limit. In this limit we found that 
the nonlinear coefficient had a positive sign. This meant 
that in the linearly stable case (62 > 0) it was possi- 
ble to be nonlinearly unstable if the initial perturbation 
was large enough. We found the critical amplitude to be 




FIG. 5. Time trace of the amplitude of density perturbations 
p (solid line) and x derivative of the density p' (dashed line) 
for b 2 /B c w 10% with a = 10" 4 . 



proportional to s/b^. 

A nonlinear numerical MHD simulation fully confirms 
the analytic result. We have used a numerical simula- 
tion of the nonlinear, full, compressible, MHD equations 
with small dissipation to verify our analytical result. We 
showed very good agreement between the simulation and 
the theory for deviations, b 2 , from B c of up to 10%. The 
numerical results show that in the short wavelength limit 
the system is nonlinearly unstable. There is some dis- 
agreement in the time evolution of the density with the 
analytical result, but this is possible since the analytic 
calculation is for an ideal system with no dissipation. We 
also discussed why a shift in B c for the linear instabil- 
ity threshold, due to dissipation, is possible at marginal 
stability and how it is harder to explain why the nonlin- 
ear result has an amplitude dependent stability. Further- 
more, the dependence is cubic so the mode grows without 
bound once it is unstable. This is even harder to explain 
as a resistive effect. 

It should be noted that the fully analytic calculation 
is facilitated by using a very simple form (a constant) for 
the transverse stabilizing magnetic field. So, while the 
conclusions of this paper seem to be on solid ground, the 
application of these findings to various systems, to the 
extent that the transverse B field of this paper is very 
special, must be appropriately qualified. For example, in 
tokamaks and stellarators, the interchange mode arises 
on rational surfaces which corresponds to a slab model 
with a sheared magnetic field vanishing at x = 0. In 
the solar coronal case, line-tying is an important charac- 
teristic absent in our simple case. Nonetheless, the con- 
clusions are sufficiently dissimilar as to indicate further 
investigation. Thus, for example, a neighboring nonlin- 
ear saturated state for the interchange mode was found 
in Refs. |S] and [7] - whereas the corresponding result in 
our case, for bi < 0, indicates a robustly growing mode 
with no nonlinear saturation. Of course, the transverse 
magnetic field in these papers was a sheared field with a 
rational surface for the unstable wave mode. Attempting 



a marginal stability analysis for sheared field, similar to 
that used in the present paper, is not straightforward. 
The fact that the sheared field goes to zero as x goes to 
zero means that a new inner ordering is required, which 
makes the calculation more involved. 

Our results are more consistent with the nonlinear in- 
stability found in Ref. |S| where the authors were also in 

— 1/2 

the parameter range with k± 3> 1, A x ~ k , ' , and 
£ x <C A x . It should be noted that their analysis was 
for the 3D line-tied g mode with no transverse field at 
marginal stability. Even so, the suprising result is that 
in both cases the system takes off once it becomes non- 
linearly unstable. This occurs even when the linear term 
is stabilizing. The primary difference between the results 
is the amplitude dependence of the nonlinear term. In 
Ref. [5] the nonlinear term has a quadratic dependence, 
while our analysis yields a cubic dependence on ampli- 
tude. If we construct an effective potential, we observe 
that the result from Ref. [8] indicates a dependence on 
the sign of the perturbation at the metastable boundary, 
while our potential is symmetric in A. Another differ- 
ence is that the result in Ref. [8|was somewhat mitigated 
by Refs. M and [TO] in that the latter papers argued that 
the ordering giving nonlinear growth would break down 
at small amplitudes before the instability fully takes off. 
In our case, our numerical simulations seem to show, in 
agreement with analytic constraints, that the nonlinear 
instability growth continues without bound and the the- 
ory only fails when A ~ 0(1) (as saturation is reached). 

Our results could also be relevant to tokamak balloon- 
ing modes to the extent that these modes are stabilized 
by an "average minimum- B well" and thus always have 
some parallel wavenumber. Work is in progress to quan- 
tify this better. Finally, our results also indicate a closer 
look at interchange stability in stellarators, presumably 
in average minimum-_B stabilized systems. 

Further investigation is necessary to answer some ques- 
tions regarding the results found in this paper. The 
transient initial growth in the time traces, mentioned in 



Sec. Ill needs to be explained. The change in the growth 
rate as the system gets closer to marginal stability, with 
£>2 < 0, needs to be investigated and compared to the 
results from Ref. [31 
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APPENDIX A 



Derivation of (|33l)-([35 



In simplifying (31 ), we found the functional, J-"[0i, 0> 2 ], 
to be 



1 9 



6 Br 

Bx • V±ViV 2 + B 2 • V±Vi^ x . (Al) 



The above equation can be simplified by writing ipi and 
■02 a certain way. From ( 24 ) we can write 



02 = 02 + 4>2 (A2) 

where 

= A{tfQ 2 {x)cos{2ky). (A3) 
Writing -02 in this way, we get the following results 

d y ^2 = dyi> 2 , (A4) 



vi0 2 



where we have written 0^ = tp\ + tpf and 
V| = ^W 2 CiW 2 cos(2fcy), 



0? = -A{t)\ 1 {xf 



(A5) 



(A6) 



(A7) 



The result ( A5 ) can be derived by multiplying ( 25 ) with 
A(t) 2 cos(2fcy) and recombining the terms. Similarly, if 
we multiply (17) by A(t)cos(ky) 7 we find that 



c 

Since B • Vj_V> = for all orders, we get that 

B x • Vi^ = 2V»i(Bx • Vj_Vi) 
= 

= B 1 -V ± (V|+02), 

and therefore 

B x • V x 0f = -Bx • Vr0f 

Similarly, since 

Bi ■ V_lV2 + B 2 ■ V_li/>i = 

= B X ■ V_l(V32 + ^ 2 ) 
+ B 2 • V_l^i 



(A8) 



(A9) 



(A10) 



(AH) 
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then it follows that 

B X • V_L^ 2 + B 2 • V.L01 = -B x • V_L^ 2 

= (A12) 



We can now simplify the last two terms in ( Al ). Using 



( A2 ) we have 

Bx • v ± vi^ 2 = Bi • v ± vi(^ 2 + 

= B r Vi(-^ 2 



S3 



■Bi ■ VxV2 



(A13) 



where we used ( A5 ) , and took advantage of the fact that 



-02 has no y dependence, to remove the Laplacians. Sim- 
ilarly, we use ( A8 1 to get 



B 2 • V ± Vi^i = B 2 • V xirjpfoh) 



(A14) 



where we used (|A4[) to get the second term. 



Combining (IA13I) and (A14) we can use (AlOl and 



(A12) to further simplify the terms with a gradient op- 



erator. So finally we get 
Bx • V ± Vi^ 2 +B 2 • V ± Vi^i = - JaVcA^ 



Bl 



Po'dylpltpf ~ dylpxlpz 



(A15) 



We can also rewrite the first term of ( Al I, 
9 



Bi 



(A16) 



We can now substitute for ipi, ip 2 i "02, ipi an d using 



(161, 030), (A3), (A6), and (TA7|. As described in Sec- 



tion III C we use the operator Jdx ^i(x) Jd(cos(ky)) on 
3l] ) in order to extract the terms that have a sin(fcy) 
dependence. The other terms will be irrelevant since the 
integration will evaluate to zero if the dependence doesn't 
match. And so we find that 



'/I /,./;) )J-r,.ro] : rrA A(t) 3 \ J^'o'Cl^ 



1 1 

~ 4 S 



(^Ci(C 1 2 )" + ^'Ci(C 1 2 )'+2Po"Ci 3 

Ci(d 2 )""l (A17) 



Finally, we use the operator Jdx £i(x) on the above equa- 
tion to get 



dx(i(x) / d(cos(ky))J r [ipi,ip 2 \ = nkL p A(t) a x 

(p'oC 2 (C 2 )"} 



9 (Ptt%) 1 ■" 



B, 



4B? 



1 1 

AB~ r 



<Ci 2 (Ci 2 )'' 



(A18) 



We made use of the fact that £1 (x) decays exponentially 
at the boundarie s to combine the three terms propor- 
tional to g / B 3 in ( Al 7 1 into one term through integration 
by parts. 



To complete the derivation of (33)-(35) we still need to 



simplify the rest of the terms. It is easy to see that after 
using the annihilation operator then we get 



dx(i(x) I d(cos(ky))(~2-^rb 2 p' d y ip 1 ) 



~2nkL p A(t)^b 2 (p / ( 2 1 ). 



(A19) 



Applying the same operator, we find that 
dx(i(x) Jd(cos(ky))B c £(ip 3 ) 

= -kB c Jdyjdx Ci sin(%) (Vi0 3 
= -kB c Jdy Jdx ^C" sin(fcy)0 3 



Bl 



P'0^3) 



+ Ci(~fc 2 sm(ky))tp 3 + Ci sin(ky)^p' ip 3 



(A20) 



—kB c Jdy Jdx sixi{ky)ip 3 X. 

(C'-^Ci + ^p'oCi), 



and therefore, using (17), 



dx&(x) d(cos(ky))B c C(ip 3 ) = 0. 



(A21) 



We, once again, took advantage of the boundary condi- 
tions to perform some integration by parts to arrive at 
the above result. Lastly, the operator on the left-hand 
side of (31 ) gives 



dx(i(x) I ' d{cos(ky))(jpp p' dyil)i - p' dyip[) 



7rkL p A(t)^{p p' tf), 



(A22) 



where the second term was thrown away since it evaluates 
to zero due to the parity of the equilibrium density. 



Collecting the terms (|A18|, |Al9j , ( |A21[ ) and ((A22j) 

together, we arrive at the (32). 
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APPENDIX B 

Description of numerical simulation 

The two-dimensional numerical simulation solves the 
following equations: 

dtp + V± • (jm±) - D„Vlp = S, (Bl) 

dt(pa±) + Vi • ( P u ± u ± ) - pVl(pu ± ) = F±, (B2) 

d t {pu z ) + V± • (pu z u x ) - pV\(pu z ) = [B z ,ijj}, (B3) 

d t B z + Vi ■ (B,u ± ) - v±^lB z = [1>,u z ], (B4) 

d t ip + u ± ■ V_l^> - ryVii/j = 0, (B5) 
where [/, h] = d x fd y h - d x hd y f and 

F x = -Vx(^P+^A-Vx^1^-P9^ (B6) 



S = VxSo f e -(.*-*i)'/*r' _ e -(x-x 2 ) 2 /2a^ _ (B7) 

The system is initialized with p = 1 and £? z = 1. We 
use Tq/M = 0.3 for the temperature and <? = 0.15 for 
the gravitational acceleration. The Gaussian function 
sources have amplitude So = 4.5, width a 2 = 6.25 x 10~ 4 
and centered around x\ — 0.7 and X2 = 0.38 (where 
L x = 1). The values are chosen by trial and error to cre- 
ate a good p' (x) profile for the simulation. The relative 



strength of the dissipation terms are as follows: 
p = rj = 5 x 10~ 4 , 

D p = 10"V 

The dissipation in the density, D p , is for numerical sta- 
bility and is made orders of magnitude smaller than the 
viscosity p. As mentioned in Sec. |III| the B z resistivity, 
rj±, is made smaller than r\ and p in order to keep B y 
approximately constant. The crossfield particle diffusion 
is set by rj±. Since the time and space scales are nor- 
malized to the Alfven speed, Va z , and the box size, L x , 
the above coefficients imply a viscous magnetic Reynolds 
number of ~ 2 x 10 3 and a Lundquist number (for mag- 
netic diffusion) of ~ 2 x 10 4 . 
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